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ABSTRACT 

Photoionization modelling allows to follow the transport, the emergence, and the absorption of 
photons taking into account all important processes in nebular plasmas. Such modelling needs 
the spatial distribution of density, chemical abundances and temperature, that can be provided 
by chemo-dynamical simulations (ChDS) of dwarf galaxies. We perform multicomponent 
photoionization modelling (MPhM) of the ionized gas using 2-D ChDSs of dwarf galaxies. 
We calculate emissivity maps for important nebular emission lines. Their intensities are used 
to derive the chemical abundance of oxygen by the so-called Tg- and /? 23 -methods. Some 
disagreements are found between oxygen abundances calculated with these methods and the 
ones coming from the ChDSs. We investigate the fraction of ionizing radiation emitted in the 
star-forming region which is able to leak out the galaxy. The time- and direction-averaged 
escape fraction in our simulation is 0.35-0.4. Finally, we have calculated the total Ha lumi¬ 
nosity of our model galaxy using Kennicutt’s calibration to derive the star-formation rate. This 
value has been compared to the ’true’ rate in the ChDSs. The Ha -based star-formation rate 
agrees with the true one only at the beginning of the simulation. Minor deviations arise later 
on and are due in part to the production of high-energy photons in the warm-hot gas, in part 
to the leakage of energetic photons out of the galaxy. The effect of artificially introduced thin 
dense shells (with thicknesses smaller than the ChDSs spatial resolution) is investigated, as 
well. 
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1 INTRODUCTION 


The main information about physical processes and the physical 
state of the interstellar medium (ISM) in star-forming dwarf galax¬ 
ies (DGs) are obtained from observed emission line spectra from 
their nebular components - gaseous nebulae or nebulae surround- 
ing the compact star-form i ng (SF) region (see the tex tbooks of 
IPopita & Sutherland II 2 OO 3 I : lOsterbrock & Ferland 11200^ . 


Line intensity ratios are used to infer the physical state (den¬ 
sity, temperature, chemical composition) of ionized regions. These 
indicators are commonly du bbed diagnost i c methods and thei r use 
date back to the seventies iSearle I fl97ll : IPagel et ^ 1979h. al¬ 
though research is still ongoing 


[Pettini & Pagel 1120041: IStasihska 


in this field (see |Pilvuginll2003l : 


I 2 OO 6 I) . As it is well known, (see 


e.g. lOsterbrock & Ferland 120051) some intensity ratios (such as 


[OIII]/14363A versus /14959A and 45007A) in range typical for 
HII regions are much more sensitive to the electron temperature 
than to electron density, whereas others (e.g. [SII]/16716A//16731A 
and [OII]/13729A//13726A) are sensitive to the electron density. In 
some cases, intensity ratios are very sensitive to both Tg and tig. In 
this case, two or more diagnostics must be used at the same time to 


determine both quantities (see lGolovatv et afll 19931 : fshaw 1[T^ . 
Physical conditions inside nebulae can thus be recovered and, from 
them, the relative abundances of ions can be determined. 

Most of diagnostic methods assume constant electron temper¬ 
atures and densities as well as ionic abundances over the whole 
ionized regime, although many attempts can be found in the litera¬ 
ture to relax this hypothesis. The assumption of non-uniformity is 
necessary for instance to explain the discrepancies between elec¬ 
tron temperatures found with different diagnostic methods. In par¬ 
ticular, the temperat ure fluctuations are characterized by the so- 
called F par ameter (|Peimb ert||l967l). A different approach was 
proposed by iMathis et al. I (Il998h . who used ratios / between 
weights of emitting region! J to ch aracterize the tem perature inho¬ 
mogeneities. Based on this work, IStasinskal l l2002t) modified the 
Peimbert’s F parameter to allow for temperature inhomogeneities 
in ionized nebulae. However, in all these methods (and in other 


* / = (N 2 n 2 V 2 )!(NiniVi), where «i and n 2 are electron densities in emit¬ 
ting regions, Vi and y 2 are their volumes, and Ni, N 2 are densities of the 
emitting ions 
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similar ones), both electron temperature and density in a single 
ionization zone are assumed to be constant, whereas this assump¬ 
tion is incorrect for most real nebulae, because density, tempera¬ 
ture and chemical abundances are a ll inhomogeneously distributed 
(see e.g. Tsamis & Peguignot l2005h for the case of 30 Doradus and 
iFrever et al. I ( l2003h for models of wind-blown and radiation-driven 
HII bubbles). Inhomogeneities are much more obvious if we wish 
to investigate an entire galaxy, where different gas phases, charac¬ 
terized by very different temperatures, densities and me tallicities, 
coexist (see e.e. lWalter et al. iboOTl : I Johnson et al, Il2012h . It is thus 
necessary to perform a more detailed photoionization modelling of 
the nebular gas (NG) in DGs that takes into account the spatial vari¬ 
ation of physical properties of the gas. We can thus obtain intensity 
maps of astrophysically relevant emission lines. Using these maps 
it is possible to calculate the synthetic spectrum of the emission 
lines along different sightlines. 

Emission-line diagnostics are also crucial to infer the correct 
star-formatio n rates (SFRs) fro m galaxies, most commonly from 
Ho- emission dKennicutt Ill998h . In fact, Kennicutt’s SFR determi¬ 
nation is based on the assumption that photons responsible for the 
hydrogen ionization stem from the stellar Lyman continuum (Lyc) 
alone and neglects other contributions as e.g. the cooling radiation 
of hot bubbles and the excitation from shock-heated regions. 

In this paper we illustrate a way to combine detailed chemo- 
dynamical simulations (ChDSs) of DGs with a state-of-the-art 
multi-component photoionization model (MPhM), to study the 
emission properties of star-forming DGs and compare them with 
observations. The idea is to use the density, temperature and chem¬ 
ical abundances obtained by means of ChDSs at different locations 
in the DG and to solve (this is the task of the MPhM) the ionizing 
radiation transfer equation as well as the ionization-recombination, 
the energy balance, the radiation diffusion equation and the statis¬ 
tical equilibrium equations in different parts of the galaxy, using as 
physical parameters for each region the results of the ChDS. An 
additional crucial ingredient is the stellar spectru m, which is pro¬ 
vided by running a tailored Starburts99 simulation dUeitherer et al. I 

I1999I1 . 


The strategy is thus to run conventi onal ChDSs and to 
use the MPhM as a pos t -process (see also Ta jma et al. 1 1201 ll : 
[Paardekooper et ani2013l : ICompostella et al.l 2013h . The obvious 
advantage of our app roach when compared to radiative hydrody- 
namical simulation s dYoshida et al. II2007I : iKrumholz et al. |[2012l : 
iHirano et al, 1[20Ti) is that the required computational time is very 
short and that not many assumptions need to be made about radia¬ 
tive transfer. The serious disadvantage of our approach is that the 
radiation field does not feedback on the hydrodynamics as it should 
and the energy and pressure in the ChDSs are purely of thermal ori¬ 
gin. 

This is the first of a series of papers in which we try to ex¬ 
tract as many useful information as possible from this combined 
use of hydrodynamical simulations and photoionization models. In 
this first paper, we will mainly focus on the details of the MPhM 
and the synergy between it and the chemo-dynamical simulations. 
Moreover, we will focus on the escape probability of ionizing pho¬ 
tons from our model galaxy, a topic of considerable relevance 
in the study of the reionization of the universe, as it is believed 
that DG-sized objects at very large redshifts are the main culprits 
for the reionization of the universe (s ee e.e. IVanzella et al.1l2012l : 
iKimm & Cenll201 j ^ lWise et al. II2014I ), although their contribution 
is still under debate (see e.g. IFuiita et al. Il2003h . Finally, we will 
in detail exploit the Ha emission from galaxies and the way it has 
been used in the past to estimate the SFR in galaxies. As men¬ 


tioned above, we will revise the estimates of SFRs based on the 
Kennicutt’s Ha -SFR relation and we will also offer a p ossible ex¬ 
plana tion for the mis match of Ho- vs. UV brightnesses dUee et al. I 
l2009h as proposed bv iRelano et al. I ( 1201 2h for galaxies with SFRs 
of less than ~ 0.01 Mq yr“*. 

The paper is organized as follows: in Section 2 the DG simu¬ 
lations are summarized. In Section 3 the details of the evolutionary 
models of the SF region are calculated; in particular, the Lyc spec¬ 
trum due to the DG star formation history is derived. In Section 
4 the MPhM is thoroughly described and the results of the pho¬ 
toionization models of different DG angular sectors are analyzed. 
In particular, we analyze the behavior of intensities of several im¬ 
portant emission lines emitted along the different radial directions 
from the SF region. In Section 5 the numerical scheme for the cal¬ 
culation of the synthetic emission line spectrum (using emissivity 
maps obtained from MPhM) for any aperture position will be given. 
Moreover, the ability of different diagnostic methods for the deter¬ 
mination of physical parameters and chemical abundances in the 
NG is tested. In the same Section we consider whether Ha is a 
good indicator of the SFR in galaxies, as it is commonly assumed. 
In Section 6 we study in detail the escape of ionizing photons from 
our model galaxy. Moreover, the change in the Lyman continuum 
spectra due to the intervening hot gas is studied in detail. In Section 
7 the most important results of this study are discussed and some 
conclusions are drawn. 


2 DESCRIPTION OF THE CHEMO-DYNAMICAL 
SIMULATIONS 

The ChDSs we use in th is work are based on the simulations of 
iRecchi & Henslerl ( l2013h . hereafter RH13. This work was aimed at 
simulating dirrs with different baryonic masses and initial gaseous 
configurations, in order to study the conditions under which a 
galactic wind could develop and determine the fate of metals, 
freshly produced during a SF episode. In all RH13 simulations, the 
galaxy is assumed to be axially symmetric and rotating about the 
axis of symmetry, with a frequency that depends only on the dis¬ 
tance from the axis R and not from the vertical coordi nate z, in 
comp liance with the so-called Poincare-Wavre theorem dTassoul I 
Il978h . The numerical code solves the 2-D gas-dynamical equations 
in cylindrical coordinates, using a second-order scheme with flux 
limiters. The evolution in space and time of various chemical el¬ 
ements (H, He, C, N, O, Mg, Si and Fe) is followed by means of 
appropriate passive scalar fields. Detailed sources of energy and of 
chemical elements due to supernovae (SNe) of both Type II and 
Type la and from winds of massive and intermediate-mass stars are 
taken into account. The details of the numerical scheme and of the 
implementation of the feedback from dying stars can be found in 
dRecchi et al. Il200ll.l2004l2006h . 

The set of models of RH13, in spite of the fact that they are 
no t tailored to any spec ific dwarf galaxy (at variance with the ones 
of iRecchi et al. I d200d) for instance), have a much higher spatial 
resolution and a more accurate description of the chemical feed¬ 
back than previous simulations of our group. They are thus ideal to 
illustrate the way in which emissivities can be calculated starting 
from ChDSs. We will use, as a reference model, the model L8F of 
RH13, i.e. a model with 10* Mq of baryonic mass and a flat initial 
configuration. This model develops a nice and undistorted galac¬ 
tic wind (see fig. 4 of RH13). Since one of the aims of our paper 
is to study the emissivities from galactic wind regions and escape 
fractions of various photons, this model is ideal for our purposes. 
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Figure 2. Syntethic radiation spectra from the SF region at different ages 
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the starting ionizing spectrum, we use tailored Starburst99 (v.7.0.1) 
simulations. The input SFR matches the one adopted in various 
ChDSs (see Section|3. The adopted initial mass function (IMF) is 
the Salpeter one. The Padova AGB evolutionary tracks for Z=0.004 
as well as the Maeder wind model were used. We run the Star- 
burst99 simulations for 200 Myr, with timesteps of 1 Myr. 

The total number of the ionizing photons Qi„„ emitted by the 
SF region per unit time weakly depends on the age. The resulting 
Lyc spectra at the inner radius of the nebular region are given in 
Figure|2] 


4 THE MULTI-COMPONENT PHOTOIONIZATION 
MODELLING 


Figure 1. Gas density (upper panel) and temperature (lower panel) maps 
obtained from the L8F model (see Sect. 2 for details on the ChDS mod¬ 
els), calculated at an evolutionary time of 140 Myr. Brighter colors indicate 
larger gas densities and temperatures. 

We have also studied the emission properties of the model L8M 
of RH13. This is a model, similar to the reference model L8F, but 
with a reduced initial degree of flattening. Consequently, a galac¬ 
tic wind develops later in this model compared to model L8F and 
in a more distorted and turbulent way (see again fig. 4 of RH13). 
In both models, the SFR is assumed to be constant during the first 
500 Myr of galactic evolution and equal to 2.67 • 10“^ Mo yr“'. In 
spite of the differences, the results obtained with model L8M are 
qualitatively similar to the ones obtained with the model L8F. They 
are only more difficult to interpret due to the enhanced turbulence, 
therefore we concentrate on L8F. We do expect significant differ¬ 
ences with the roundish model L8R though. A detailed comparison 
between emissivities and spectra obtained with models L8F, L8M, 
and L8R will be presented in a forthcoming paper. 

A typical map of density and temperature resulting from the 
ChDS of the model L8F is shown in Figure[T] 


3 SPECTRA OF THE IONIZING RADIATION FROM 
STARS 

The Lyman continuum spectrum originates from the stellar popu¬ 
lation and then is modified by the intervening gas. In order to set 


In this section we present the algorithm to perform the multi- 
component photoionization modelling of NG in DGs. The results of 
such modelling based on ChDSs, such as the emissivity maps for 
important emission lines, the transformation of the ionizing Lyc- 
spectrum and the predicted emission line spectra obtained along 
various radial directions will be considered as well. 


4.1 The MPhM wrapper 

As basis for the photoioni zation modelli ng we used G. Ferland’s 
code CLOUDY 08.00 (see lFerlandll2008h . 

Looking again at the results of ChDSs (see Fig. [T]l, the first 
important conclusion we can draw is that the ionizing gas can be 
divided into two main parts. The first one is the cavity carved by 
supernova explosions and stellar winds around the SF region (su¬ 
perwind region - SWR). This component contains very rarefied and 
hot gas. Also the abundance of heavy elements is very high. Since 
the MPhM does not treat shock-heated gas, the temperature in this 
region can only be determined by the ChDS. The second compo¬ 
nent is the supershell swept up by the propagating shock (hereafter 
called ’wall’). It is characterized by much higher gas densities, as 
well as by heavy-element abundances typical for the galactic ISM. 
The ’wall’ corresponds thus to the outer edge of the SWR. It is 
important to point out that, in spite of the larger densities and the 
relatively low temperatures, a significant fraction of the wall is ion¬ 
ized because the leakage of ionizing photons from the inner region 
of the galaxy is significant (see below). Beyond the wall, the gas is 
unperturbed and it is heated mainly by photoionization. Therefore, 
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the temperature in this component can be obtained as a solution of 
the photoionization thermal energy-balance equation. 

The second important conclusion from the analysis of ChDS 
results is that the maps of density and temperature are very com¬ 
plex. In spite of the axial symmetry, the whole galaxy can not be 
approximated by a simple geometrical figure. Therefore, it is im¬ 
possible to use any of the geometrical shapes offered as options in 
CLOUDY. Obviously, a more complex strategy to simulate the gas 
emission is required, in order to take into account the very large 
variations of density, temperature and chemical composition across 
the galaxy (see also th e Introduction for alternative strategies). 

iMorissetl ( l2006h (see also references therein) developed the 
so-called pseudo-3D extention Cloudy JiD of CLOUDY, in order 
to model the emission of nebular objects in 3D geometry, using 
many ID photoionization models. This is an appropriate strategy 
to model objects like diffuse ionized gas surrounding star-formation 
regions. Nevertheless, we decided not to use this code and to write 
our own specialised driver MPhM (multicomponent photoioniza¬ 
tion modelling) to calculate the multicomponent photoionization 
models for the following reasons: 1) We are not using the standard 
CLOUDY release but a customized version. In particular, it was 
necessary to update the core of CLOUDY to distinguish between 
SWR and ordinary ne bular gas. 2) To g enerate brightness maps end 
emission line profiles. lMorisset I d^Om developed IDL-procedures, 
but IDL is not a free software and therefore it is not available to all 
astrophysicists who are interested in such models. Therefore, we 
developed our own code DiffRay (described below) to calculate the 
emission line spectra in synthetic apertures with user-defined sizes 
and positions. We plan to make the MPhM driver as well as the code 
DiffRay freely available to all the astrophysical community. 3) We 
are continuously upgrading MPhM (and also the core of CLOUDY) 
to better adapt it to the study of the nebular emission in galaxies. 

In the MPhM approach, the whole computational box of the 
ChDS is divided into angle sectors, namely solid angles drawn from 
the origin of coordinates, which is the location of the central star 
cluster (see Figure^. As already mentioned, ChDSs are calculated 
in 2D geometry assuming cylindrical symmetry. This implies in 
particular that only one quadrant in the R - z plane must be calcu¬ 
lated. We divide this quadrant into 20 angle sectors (see Figure^. 
In what follows, the sectors will be numbered from 1 to 20, with 
Sector 1 being the one adjacent to the galactic plane and Sector 20 
the one adjacent to the z-axis (the axis of symmetry). Each sector 
is in turn uniformly divided into 200 radial components. The radial 
extent of each component is 11.5 pc. The gas at radial distances 
larger than 2.3 kpc has been considered as part of the intergalactic 
medium. 

A simple coordinate transformation and an interpolation allow 
us then to map the ChDS results into the MPhM frame. In particu¬ 
lar, in each of the sectors and components, the following input data 
have been derived and stored 

• the hydrogen density; 

• the electron temperature (only in the SWR, i.e. in the inner¬ 
most part of sectors); 

• the chemical abundances of 9 chemical elements (H, He, C, 
N, O, Mg, S, Si, Fe). 

It is worth stressing that CLOUDY itself divides the computa¬ 
tional volume into zones which are usually smaller than our above 
defined components. Therefore, CLOUDY performs a further in¬ 
terpolation on the input data. 

The edge between hot SWR and wall is defined by the locus 
of points across each sector where the density gradient is highest. 
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Figure 3. Scheme of the division of the computational box into sectors and 
components. Sectors are numbered from 1 to 20, with Sector 1 being the 
one adjacent to the galactic plane (i.e. adjacent to the R-axis) and Sector 20 
the one adjacent to the axis of symmetry (i.e. adjacent to the z-axis). 

In most of the wall (and also in external regions), the temperature 
derived from the ChDS is below ~ 4000 K. For these components, 
we rely on the temperature derived by CLOUDY, as the heating 
effect due to radiation can be very significant. A detailed analysis 
of such situation is presented in the following subsection. 

In our models, two kinds of ionizing radiation are taken into 
account. The first is the stellar emission from the SF region. In 
Fig-Elthe synthetic ionizing radiation spectra are shown, depending 
on time. The second kind of ionizing radiation is the diffuse one ap¬ 
pearing in the NG, mainly due to the radiative recombination onto 
the ground level of the ions H^, He^, He^^, onto the second level 
of He^^ as well as by Lya transitions of Hel and Hell. The diffuse 
ionizing radiation is calculated using the ”outward-only” approx¬ 
imation, which is the default in CLOUDY. According to this ap¬ 
proximation, the diffuse radiation flux from each component, (tak¬ 
ing into account all important sources of opacity), is added to the 
incident continuum flux from the central SF region. The radiation 
from each component is allowed to propagate only outwards. Thus, 
each angular sector in the present MPhM approach is ’’isolated” 
from the ionizing radiation from the others. This is clearly a sim¬ 
plification, made in order to reduce the computational costs. More 
advanced solutions are currently being developed and tested. 

We perform MPhM calculations in the time interval from 10 
to 200 Myr with timesteps of 10 Myr. At each time step, we take 
a different ChDS output and a different stellar spectrum, in order 
to keep pace of the changes occurred in the gaseous and stellar 
component of the galaxy. 

4.2 Photoionization modelling of sectors 

As mentioned in Section [2] we will mostly focus on MPhM cal¬ 
culations based on the ChDS L8F. In Fig. |4] the variation of the 
electron density and temperature across some angular sectors are 
drawn, as a function of radial distance from the center. Different 
lines indicate different evolutionary times. As a reminder. Sector 1 
runs across the galactic plane, therefore the hole of low density and 
high temperature (the SWR in our nomenclature) is less extended. 
Conversely, Sector 19 lies close to the axis of symmetry, i.e. along 
the direction of steepest pressure gradient. Along this direction a 
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Figure 4. Variation of the electron density and temperature across some angular sectors, as a function of radial distance from the center. Different lines indicate 
different evolutionary times. The calculation refers to the ChDS model L8F (see Sect.|2]for details). 
























































































6 Melekh et al. 


galactic wind develops faster and thus the SWR is more extended 
(see also Fig.[T). 

In Fig.Uone can clearly see the edge between the SWR, char¬ 
acterized by high temperatures and low densities, and regions of gas 
still at lower temperatures and at higher densities. In some cases 
however some oscillations in the temperature profiles can be dis¬ 
cerned, with colder regions preceding hotter ones. This is due to 
the turbulent nature of the galactic outflows and to the adiabatic 
cooling of the outflowing gas. This can be clearly seen in Fig.[T] 
A ridge of shock-heated gas stretches from (R,z) (0, 2) kpc to 

- (0.8, 1.5) kpc. This gas has a temperature higher than the back¬ 
ground one. Behind it, there is a strip of adiabatically cooled gas. 
Moreover, a turbulence-generated eddy starts developing at (R,z) ^ 
(0.5, 1.4) kpc. It is thus possible that some sectors, at some evolu¬ 
tionary times, intercept regions at different temperatures, generat¬ 
ing thus the temperature oscillations visible in Fig. |4] This effect 
can be particularly significant from Sector 10 onwards. Because 
of these inhomogeneities in the temperature distribution, we allow 
CLOUDY to calculate the electron temperature any time the tem¬ 
perature derived from the ChDS drops below 4000 K. It can also 
be seen that for sector 1, containing the densest and thickest ’wall’ 
element, the outer ionization front occurs at distances from 0.8 to 
1.3 kpc, depending on age. This is the outer ionization front due to 
the absorption of ionizing photons in the gas. This means that in 
radial directions of low-number sectors the gas opacity is large and 
a very small fraction of the ionizing radiation propagates beyond 
the wall. 

The MPhM allows to calculate the emissivity maps for impor¬ 
tant nebular emission lines, that can be used to calculate the corre¬ 
sponding radiative fluxes along any direction. This is necessary in 
order to obtain theoretical spectra to compare with observed ones. 

In order to investigate where in the nebula the recombination 
as well as collisionally excited (forbidden) emission lines arise, we 
have plotted the radial distributions of line emissivities (see Fig- 
ures[5]-[6}. We have chosen to plot only four representative sectors 
and not all 20 in order to make the plots clearer. In particular we 
have chosen the plane sector (sector 1), two central sectors (sec¬ 
tors 8 and 13) and a sector close to the axis (sector 19). We have 
not chosen sector 20 because axysimmetric hydrodynamical simu¬ 
lations (such as our simulations) unavoidably show some numerical 
noise along the symmetry axis. In the left panels of Figure |5] the 
radial distribution of the Ha emissivity (as usual, separated into an¬ 
gular sectors) is shown. This line emissivity behaves quite similarly 
to the electron density (see Figure|4l( as expected, because this re¬ 
combination line has a strong dependence on electron density and 
a very weak dependence on the electron temperature. For the first 
seven sectors the emissivities drop at the ionization front, whereas 
for high-number sectors the external regions of the galaxy are emit¬ 
ting significantly in Ha . This is because the optical thickness of the 
gas is lower and the radiation field (mainly due to Lyc photons from 
the central SF region) keeps the temperature until R ~ 1 - 2.5 kpc 
relatively high, with T,, ranging between 4000 and 32000 K. The 
optical thickness for Lyc photons is not high, because the interven¬ 
ing gas has on average relatively low densities (of the order of 0.1 
cm“’). This figure is very low in comparison with typical densities 
in conventional HII regions, where iih = 1 - 300 cnT^. Notice that 
in spite of the good spatial resolution of our hydrodynamical simu¬ 
lations (4 pc), this is still not enough to resolve HII regions properly 
(let alone to resolve compact and ultra-compact HII regions) and 
this is clearly a limitation of our model, as the leakage of Lyc pho¬ 
tons in our results is perhaps higher than it should be (see Sect. [U. 
High-resolution multi-phase simulations are currently being tested. 


in order to remedy to this situation. In this work, for illustrative 
purposes, we have decided to make a simpler (and less physically 
justifiable) hypothesis, namely that a clumpy phase is present in 
the galaxy at spatial scales that can not be solved by our hydrody¬ 
namical code. Therefore, we add this clumpy phase only during the 
MPhM post-processing (see next Section for more detail). 

The behaviour of \{p and HeI/14471 as well as other Hel lines 
emissivities are also quite similar to the Ha . The emissivity of 
HeII/14686 shows the peaks in ’wall’ sectors, corresponding to 
maxima in the electron density distribution. Again, for sector 1, 
this emissivity drops to almost zero at the ionization front. This 
is because the ionization potential of He+ is 4Ryd. The ’wall’ ab¬ 
sorbs the photons beyond 4Ryd. Therefore, the fraction of the ion¬ 
izing photons beyond the ’wall’ is very small. Thus, at the inner 
edge of the ’wall’ one has frequently He++/He>He+/He, but the 
He++/He fraction drops faster than the He+/He one. The emissivity 
of HeII/14686 line is only caused by recombination of He^^. In the 
wind region (SWR) He++/He»He+/He, because the high tempera¬ 
ture here prevents the recombination so that there does no signifi¬ 
cant emissivity in HeII/14686 line emerge. 

In order to illustrate the behaviour of forbidden line emis¬ 
sivities we show in Fig. the radial distribution of [OIII]/15007. 
This line, like other forbidden lines, is sensitive to the electron 
temperature. The highest temperature values occur inside the su¬ 
perwind cavity region. Nevertheless, the [OIII]/15007 emissivity is 
very weak, because at such high temperatures the abundance of 

ions is very low (see next Section). As expected, for sector 
1 , characterized by the presence of an outer ionization front, the 
[OIII]/15007 emissivity drops to almost zero at the ionization front. 


5 SYNTHETIC SPECTRA, THEIR ANALYSIS, AND 
COMPARISON WITH OBSERVATIONS 

Observed emission line spectra from spatially distributed NG sur¬ 
rounding the SF re gion are usually obtained from d ifferent aperture 
positions (see e.g. iKobulnickv & Skillmanlll997L KS97). Across 
these apertures, abundances of various chemical elements can be 
derived by different diagnostic methods. The aim of this section is 
to compare the chemical abundances directly obtained from ChDSs 
with the ones obtained after applying diagnostic methods to the 
MPhM output files. In this way we try to assess the validity of var¬ 
ious diagnostic methods, as well as the correctness of our meth¬ 
ods and procedures. In fact, as illustrated in the previous Sections, 
our combined ChDS-rMPhM simulations allowed us to produce 3D 
maps of emission line and continuum emissivities, which we illus¬ 
trated by means of radial variations across various angular sectors. 
On the other hand, by means of ChDSs alone we have the spatial 
distribution of the abundance of nine chemical elements (H, He, C, 
N, O, Mg, Si, Fe), as a function of time. We thus have all required 
information to test the diagnostic methods. 

In this Section we shortly describe the procedure to calculate 
the synthetic emission line spectra from different aperture positions 
using emissivity maps. We firstly define non-central (i.e. along 
lines-of-sight which do not pass through the center of the DG)) and 
central aperture positions (see Fig.|7j. Synthetic spectra are calcu¬ 
lated by integrating the intensity maps across each aperture (see 
next subsection for more details). Synthetic spectra are in turn used 
to determine the oxygen abundances using two well-known and 
wid ely used diagnostic methods, nam ely the two-zones -method 
(see lPagel et al.Tl992l : lGamettll992h and the R 23 -method. The lat¬ 
ter is based on calibrating relationships O/H-R23 from iMcCaughI 



























8 


Melekh et al. 


Sector 1: 



Sector 13: 



Sector 8: 



Sector 19: 



Figure 6. Evolution of the radial distribution of /i5007[0 Ill] emissivities for selected sectors and ages. 


lll99ih . Then we present the comparison of the obtained oxygen 
abundances with the ones from the ChDS for the corresponding 
aperture. 


5.1 Synthetic spectrum calculation 

To calculate the emission line spectrum using emissivity maps we 
developed the 3D code DiffRaY. This code integrates the fluxes in 
emission lines over the solid angle that is defined by the aperture 
position relative to the observer position. 

We consider each aperture not as an infinitesimally thin nee¬ 
dle, but as a geometrical figure having a section of 49 pc x 49 pc. 
This section area is further subdivided in u x n smaller volumes. 
The subdivision process is iterative and n represents the stage of 
the iteration, n is equal to four in the first iteration. At each stage 
of iteration, the radiative fluxes are integrated across each subdivi¬ 
sion of the aperture. Among different possible numerical integra¬ 
tion techniques fsee lPress et al7iri992h . we have chosen the trape¬ 
zoid method;_asJH^_the_optimal choice for this kind of problems 
(see lBuhaienko & Melekh 1 20 13h . The code repeats this procedure, 
increasing u by 2 at each step, and iterates until convergence, i.e. 
until the differences between line intensities at two different stages 
of iteration are less than 2 %. 

We tested this code by comparing the radial line fluxes ob¬ 
tained from this 3D integration program with the ones obtained by 
CLOUDY. In Table[T]the ratios between our results and CLOUDY’s 


ones are given. It can be seen that these ratios are always very close 
to unity, with deviations of ~ 1.7% at most. 

5.2 Synthetic spectra from different aperture positions 

In order to test the ab ility of the two-zo nes T,,-method jPagel et al. I 
1 1 992 I : lo^nett II 1 992h and 7^23-method jPagel et al, Ill979ll to deter- 
mine the oxygen abundances in galaxies, we calculated synthetic 
emission line spectra using MPhM emissivity maps, as described 
above. 

The angle a = 30“ (see Fig|7} characterizing the inclination of 
the galaxy relative to observer was adopted. 

5.3 Comparison of modelled and observed emission line 
spectra 

The aim of this subsection is to compare the calculated relative 
line intensities with the observed ones. For this comparison with 
observations we use the sample of selected observed spectra from 
iKehrig et al. I l l2004h . From this sample, we selected galaxies with 
Ho'/FljS between 2.6 and 3.2, because our models results are in this 
range (see Figure[8l(. The resulting observational ranges for the in¬ 
tensities of important emission lines are given in column 2 of Table 
Of course, we do not expect a perfect match between theoretical 
expectations and observations, because this paper is not intended 
for this purpose, but model emission line intensities should be the 
of the same order as observed ones. 
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a 



side view face on 


Figure 7. The definition of aperture cuts. The central aperture is labeled by zero, a is the angle between the symmetry axis and the plane normal to the observer 
sight line. The left panel is the side view of the object. The right panel is the front view. Here, the size of the aperture can be appreciated. 


Table 1. The ratios between radial fluxes, calculated by our 3D integration 
program DiffRaY and the ones obtained by CLOUDY for the sectors 1, 8, 
13, and 19. 


Line 

Sector 1 

Sector 8 

Sector 13 

Sector 19 

HJ3 

0.985 

1.002 

1.002 

0.985 

Ha 

0.985 

0.994 

0.997 

0.995 

Lya 21216 

0.984 

1.011 

0.991 

0.983 

[OII]43729 

0.998 

1.003 

0.999 

0.990 

[OII]43726 

0.999 

0.992 

0.998 

0.998 

[OIII]/i5007 

0.984 

0.990 

0.996 

1.003 

[OIII]/i4363 

1.000 

0.994 

1.000 

0.979 

Hel 44471 

0.998 

0.991 

0.989 

0.993 

Hel 45876 

0.994 

0.991 

0.995 

1.004 

Hel 46678 

0.995 

0.996 

1.006 

1.002 

Hel 47065 

0.987 

0.991 

0.993 

0.994 

[SII]46716 

0.988 

1.002 

1.008 

0.984 

[SII]46731 

1.004 

0.992 

0.993 

1.000 

[NII]46548 

0.989 

0.995 

1.004 

1.001 

[NII]46584 

0.989 

0.996 

1.004 

1.000 


Table 2. Comparison of relative line intensit ies observed at galaxies, se¬ 
lected from the sample of iKehrig et al. I i2Q04l) with derived ones of models 
Ml and M2 (see text). 


Relative intensity 

Observed 

range 

Ml 

M2 

[OIII]45007/Hj3 

0.58. 

. 7.84 

2.08 

3.47 

[OIII]44363/Hj3 

0.03 . 

. 0.15 

0.03 

0.03 

[OII]43727/H/3* 

0.40. 

. 7.19 

0.18 

1.95 

[OIII]45007/[OII]43727 

0.17 

..4.8 

11.6 

1.78 

[SII]46716/H8 

0.05 . 

.0.66 

0.02 

0.23 

[SII]46731/H8 

0.04. 

.0.76 

0.02 

0.16 

Ha/HJ3 

2.63 . 

. 3.18 

2.95 

2.93 

[NI1]46584/H« 

0.02 . 

.0.29 

0.001 

0.025 

[SII]46716/Hq' 

0.02 . 

. 0.21 

0.007 

0.08 

[SII]46731/Hq- 

0.02 . 

.0.25 

0.007 

0.06 

[SII](46716 + 46731)/Hq- 

0.03 . 

.0.46 

0.01 

0.13 


"[OII]/13727/H/3 = [0111(23726 + 23729)/H^ 
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Figure 8. The evolution of Ha/H^g for six synthetic aperture positions (AP). 


In Table 2 and in what follows, we label with Ml the standard 
model, i.e. the analysis done starting from the results of the ChDS. 
Model M2 refers to a modification of the ChDS results, in which a 
very thin shell of enhanced density is artificially added (see below 
for more detail). 

In Figuresl^the time dependences of the relative line intensi¬ 
ties (relative to Hj3) of important emission lines, obtained for var¬ 
ious synthetic aperture positions (see Fig. |7) are shown. It can be 
seen that our models predict [OIII]/15007 (see Fig.|9^)) line intensi¬ 
ties higher than the observed minimum or very close to the one for 
most of APs. Only in the case of AP 4 intnesity of this line is higher 
than the observed minimum after 120 Myr. The relative intensities 
of [OIII]d4363, [OII]T3727 and [SII]d6716 (see Fig.|9l3-d) ) for 
most of ages and APs are not able to reproduce the observed mini¬ 
mum. 

To further simplify the comparison of the model results with 
observed ones, we select the relative intensities obtained as a 
function of radial distance of sector 1 at 140 Myr. The result¬ 
ing intensities are given in column Ml of Table [2] This angu¬ 
lar sector has the same problems shown by the aperture positions 
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0 20 40 60 80 100 120 140 

Age, Myr 


c) 



b) 




Figure 9. Some important line intensities, obtained for six synthetic aperture positions (AP), as a function of time: a) [OIII]/15007 - one of the main coolers of 
nebular gas, also used in diagnostics for oxygen abundance determination; b)[OIII]/14363- important line for electron temperature determination; c) [OII]/13727 
- it is a sum of [011]/13726 and [011]/13729 line intensities, that are important to determine the oxygen abundance; d) [SII]/i6716 — the ratio of this line to 
[Sll]/i6731 is used in diagnostic for electron density determination. The intensitie s are normalized to /ha- To compare with observations the minimal value, 
obtained from observed emission line spectra of galaxies, selected from sample bv iKehrig et al. I l2004ll are shown by solid lines (Min.obs.) in each panel. 


above, i.e. it severely underpredicts the [OIIJ/13727 and [SII]/16716 
emission line intensities. We also compare the intensities of 
Ha/Hye, of [NII]T6584 /Hq', and of the sum of [SII]46716/Ha and 
[SII]/i6731 /Ha with observed ones. It can be seen that only the val¬ 
ues of [OIII]T5007/H/?, Ha/Hj8 and [OIII]T4363/I^ are within the 
observed ranges. The other line ratios are underestimated in model 
Ml. 

The difficulty in reproducing the correct [Oil] and [SII] inten¬ 
sities boils down in the end to the relatively low gas density across 
most of the computational volume; a limitation of the model we 
have already outlined in the previous Section. Due to this low den¬ 
sity, the temperature beyond the ’wall’ in most cases is relatively 
high. This prevents the formation of a significant fraction of 
ions. In this part of the galaxy, the oxygen is in fact mostly doubly 
ionized and this explains the reasonable agreement between [OIII] 
intensities and observations. To illustrate this point, the various ion¬ 
ization stages of oxygen across the wall at t=140 Myr are shown in 
Fig. [To] upper-left panel. We have chosen the angular sector 1 (i.e. 
the one adjacent to the plane of the galaxy, see previous Section), 
between 0.5 and 0.75 kpc because the radial distribution of the var¬ 
ious ionization stages can be better appreciated. 

In order to address the problem of the lack of very dense re¬ 
gions in our computational volume, we decided to introduce one 
dense shell of enhanced density (TDS - thin dense shell). This 


TDS represents a thin density enhancement, due to a shock, im¬ 
possible to resolve with the currently adopted numerical resolution. 
The shell extends for 10 pc (see Fig. lilt and reaches a peak den¬ 
sity of iiH = IcirT^. This peak density (not the entire shell; just 
the peak) extends for 1 pc. The Ml density at the peak position is 
instead about 0.2 cm“^. With such a background density, the TDS 
peak density would thus represent a density enhancement of ~ 35 
in a shock. For a strong isothermal shock, this would cor respond 
to a M ach number of 6 (for instance see the textbook of lPrainel 
1201 ih . perfectly reasonable under the studied conditions. 

To illustrate the effect of the TDS presence we have chosen 
again Sector 1, because it is the sector adjacent to the plane of the 
galaxy, where the presence of a TDS is more likely. The TDS in¬ 
herits the chemical composition of the underlying material. We cal¬ 
culated models of Sector 1 that differ by the position of this TDS: 
at the wall; midway between the wall and the outer radius of the 
considered computational box (7?„„, 2.3kpc); and very close to 

Rout- Only in the case of a TDS positioned in the wall we obtained 
relative intensities in agreement with observations. The results of 
this model, labelled M2, are shown in Table 2. Also, in Figure [TOl 
(upper-right panel) the various ionization stages of oxygen for M2 
are shown. The TDS produces a significant increase of the 0^/0 
ratio at a distance larger than 0.55 kpc. Now, the radial profile of 
of in model M2 is more complicated (there are two maxima). 
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a) Model M1: L8F, Sector 1, 140Myr 


b) Model M2: L8F, Sector 1, 140Myr 


Q 

O 



0.5 0.55 0.6 0.65 0.7 0.75 

Radius, kpc 
c) Density distributions 



0.5 0.65 0.6 0.65 0.7 0.75 

Radius, kpc 


d) Temperature distributions 



Figure 10. Radial distribution of the oxygen ionization fractions, hydrogen densities and electron temperatures in Sector 1 at 140 Myr in models Ml and M2. 
Model M2 is characterized by the addition of a thin shell of enhanced density in ’wall’ with a peak density of IcmT^ (for comparison, the standard density 
obtained from ChDS at this position is about 0.2cm“^). 


This is enough to significantly increase the [O II] and [S II] emis¬ 
sion. The disagreement with the observations disappears for M2 
(see again Table 2). The values of /utr/Zus, /r 6584 [NII]//Ht, as well 
as (/u 67 i 6 [SII] + / 4673 i[SII])//Ha are within the observed range too. 

The combination of / 45 oo 7 [OIII]//h^ and / 46584 [NII]//Ha al- 
lows us to place ou r models in the so-called BPT diagram 
([Baldwin et al. I IT98lh . the famous seagull-shaped diagram able 
to distinguish between HII regions due to star formation and 
AGNs. For model M2 (log /r 6584 [NII]//H. , log /,5n n7[OIII1//H«) = 
(-1.14 ,0.54). Comparing these values with fig. 5 of iBaldwin et al. I 
( Il98ll) . it is clear that model M2 is well within the HII regions. Our 
results obey pretty precisely eq. 8 of the Baldwin et al. paper. 

In the SWR, the oxygen is present in and even higher 

ionization stages, because of the high temperatures. The density 
and temperature distributions for models MI and M2 as well as 
the temperature level of the ionization front are shown in Figure 
[T^ (panels c and d). Because of the TDS, the temperature drops 
significantly below the ionization front, at distances of the order of 
0.7-0.9 kpc. At larger distances, the temperature increases again, 
because of the reduced density. By means of the TDS, we are thus 
able to create a layer of neutral gas between two ionized regions. 

In the Ml model, the oxygen at large radii exists mainly in 
the ionization stage, therefore the abundance is too low 
to reproduce the observed [Oil] line intensities. The presence of 
the TDS (model M2) leads to the appearance of two regions of 
enhanced by fractional abundance: the first one is a thin region 


at the inner edge of the wall; the second one extends up to 0.73 
kpc. These shells are enough to reproduce the observed [Oil] 
emission in observed range (see the results of model M2 in Table 
2 ). 

Although the results of model M2 are encouraging, it is impor¬ 
tant to point out that our results are quite sensitive to the TDS posi¬ 
tion and also to the density profile within the TDS. It must be also 
emphasized that, in our models, diffuse ionizing radiation (emitted 
by ionized diffuse gas) is calculated in outward only approximation. 
Therefore, as already outlined, ionizing photons from neighbor sec¬ 
tors cannot penetrate. Therefore, the propagation of the ionizing 
radiation extends up to the outer ionization front, that occurs at » 
1 kpc. Such a restriction for the propagation of ionizing photons 
does not exists in real galaxies. We are working on an extension of 
the MPhM code to deal with this kind of non-radial propagation of 
ionizing photons. Such extension is especially necessary to include 
compact, non-axisymmetric clumps in the model, and to consider 
spatially extended and non-spherical sources of ionizing radiation. 


5.4 Diagnostics of synthetic emission line spectra 

Our study aims also at comparing the derived oxygen abundances 
with the ChDS ones, by means of two popular diagnostics, the 7’,,- 

method and the R 23 -niethod. __ 

According to 7’^-method (see jPagel et ^ 1 19921 : [Garnett I 
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Figure 11. Modification of density distribution in model M2 in comparison 
with standard model Ml calculated at 140Myr. 


1 19921) . for the determination of the oxygen abundances, we also 
need the relative intensities of [Olll](/i4363,/14959,/l5007), and 
[01I](/13726,/i3729), obtained from the synthetic emission line 
spectra and calculated above for different aperture positions. 

However, frequently the emission line 44363 [OIII] is absent 
in real observed spectra, because of low intensity. In such cases 
the so-called calibration methods are used. In such methods the 
relationships between (12 + log O/H) and the parameter /?23 = 
(f[oii] 33726 +r 3729 + f[oiii] 24959 + 25007 )/fH^ are determined on the basis 
of appropriate calib rations. We select for our tests the calibration 
0f lMcCaugh1 ( ll99lh . 

In Table 3 we compare the oxygen abundances, the 
luminosity-weighted (L-W) and mass-weighted (M-W) with the 
ones calculated for the six selected synthetic apertures, using the 
chosen diagnostic methods, for three different ages. The behaviour 
of O/H abundances along sight lines becomes more complicate 
with age. The Tf-method reproduce the L-W values within an er¬ 
ror (A) of only 0.5 dex for three APs at 10 Myr, for one AP at 
70 Myr and for two APs at 140 Myr. Larger deviations from L-W 
results can be explained by the inhomogeneous distribution of 
and ions. Also, the 7^23-method does not reproduce the oxygen 
abundance in AP 5. In most cases, the /? 23 -method overestimates 
the oxygen abundances with deviations of more than 0.5 dex, but 
for some apertures, the L-W abundances are very close to the ones 
calculated with the 7?23-method. 

These results show that the considered diagnostic methods 
cannot be used to precisely determine the oxygen abundance in 
DGs due to both the inhomogeneity of the element abundance and 
the very low density (<3cm“^ in some cases). The results show also 
that the model limitations, already pointed out above, might render 
our model galaxies quite different from real galaxies. In particular, 
as already pointed out, due to the low density, the optical thickness 
for the emission lines along most of sight lines are low (probably 
lower than in real galaxies). 

Also in Table 3 (last two rows) we compare L-W and M-W 
oxygen abundances with the ones determined using both T,,- and 
7?23-methods for models Ml and M2 at 140 Myr. It can be seen that 
for Ml the difference between the T,, and the 7?23 method amounts 
to 0.4 dex. The differences between the L-W value and the results 
of Tn and 7^23 methods are slightly larger (a bit more than 0.5 dex). 
Model M2 shows instead a good agreement of the oxygen abun¬ 
dances determined by diagnostic methods with both L-W and M-W 


Table 3. (12 + log O/H) values obtained for the ChDS model weighted 
by luminosity (L-W) and mass (M-W), respectively, and compared to abun¬ 
dances derived by means of the Tg and the 7?23 diagnostic methods. The 
two lowermost lines give the (12 + log O/H) abundances for both analytical 
model distributions Ml and M2 at 140 Myr, respectively. 


Age, 

AP 

Weighted, from ChDS 

Diagnostic methods 

Myr 


L-W 

M-W 

Tg 

7?23(McCaugh) 

10 

Central 

7.90 

7.80 

6.90 

6.48 

10 

1 

7.14 

7.23 

7.16 

6.80 

10 

2 

7.70 

7.58 

121 

7.01 

10 

3 

7.89 

7.74 

6.96 

8.98 

10 

4 

5.78 

5.78 

5.84 

6.91 

10 

5 

5.78 

5.78 

7.01 

6.67 

70 

Central 

8.52 

8.38 

7.38 

6.72 

70 

1 

8.29 

8.29 

6.84 

6.61 

70 

2 

9.13 

9.09 

6.96 

9.06 

70 

3 

8.51 

8.37 

7.48 

9.01 

70 

4 

5.77 

5.78 

5.84 

6.92 

70 

5 

5.78 

5.78 

6.95 

6.73 

140 

Central 

8.93 

8.98 

8.05 

8.91 

140 

1 

7.79 

8.22 

7.42 

6.97 

140 

2 

8.51 

8.53 

7.44 

8.99 

140 

3 

9.00 

9.03 

7.91 

8.94 

140 

4 

8.44 

8.43 

8.35 

8.82 

140 

5 

5.78 

5.78 

7.38 

6.88 

140 

Ml 

8.07 

7.78 

7.51 

7.11 

140 

M2 

8.21 

7.79 

8.10 

7.96 


determinations. In this case, the Tg and 7?23 methods can be used for 
the determination of the oxygen abundance. 

Thus, we can conclude, that Tg- and 7723-methods have limita¬ 
tions in two circumstances: (f) in the low-density limit in nebular 
enviroments with complex element abundance distributions along 
sight line; (ii) when an additional radiation field contributes to the 
stellar radiation, as e.g. the cooling radiation of a hot plasma like 
the superwind bubble. It is clear that a more detailed and punc¬ 
tual analysis of the diagnostic methods can be done only once our 
models are able to reproduce more faithfully the real physical char¬ 
acteristics of galaxies. Notwithstanding these limitations, we be¬ 
lieve that the preformed test is reasonable - we try to analyze the 
galaxy with the “eyes of an observer” and try to recover the galactic 
properties from observational methods. At some level what we “get 
out” should match what we “put in”. Moreover, since we use two 
different diagnostic methods, we expect at least some agreement 
between them. The failure to match these expectations entitles us 
to conclude that diagnostic lines have limitations in low-density, 
chemically complex regions. 


5.5 The empirical determination of the star formation rate in 
galaxies 

We conclude this section with an estimate of the reliability of em¬ 
pirical indicators for the star formation rates in external galaxies. 
One of the most commonly used proxy of SF in galaxies is the Ho- 
emission. In particular, the calibration of iKennicutt I (Il998h 

, UHa) 

SFR(Meyr-‘) = (1) 

1.26 X 10^^ ergs s ^ 
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Figure 12. The ratio of the star-formation rate determined from L(HQr) (see 
Eq.B and the ’true’ rate, i.e. the rate which is used as input from our ChDS 
(2.67 ■ 10-2 Mo yr-‘) of model Ml. 


(based on the Salpeter IMF) is very widely used. Since we are able 
to calculate the Her luminosity by means of the MPhM code, we are 
in the condition to compare the SFR based on L(Hcr) with the ’true’ 
SFR, i.e. the rate which has been used as input for our ChDS (equal 
to 2.67 • 10-2 yj.-i^ ggg Sect.O. This comparison is shown in 
Fig-im 

We can see from this figure that, at the beginning, there is a 
perfect match between the Her -based and the ’true’ SFR. This re¬ 
assures us about the correctness of our procedures and the accuracy 
of the Starburst99 code. At later times, however, the Her luminosity 
increases due to the contribution of the warm gas. The Kennicutt 
SFR estimates becomes thus inaccurate. It overestimates the true 
SFR by up to 60 % at t=30 Myr. After 160 Myr, Lyc photons begin 
to leak out of the galaxy (see Sect. [Hi, thus the global Her luminos¬ 
ity decreases and the relation Eq.[T] starts underestimating the true 
SF rate. An accurate study of the status of the gas in a galaxy is 
thus a fundamental prerequisite in order to assess the reliability of 
the SF rate estimate based on Htr emission. Given the model limita¬ 
tions, an agreement between model results and expectations within 
a factor of two should be considered as fair. Still, we believe it is 
important to quantify these inconsistencies and to outline which 
physical processes can lead to inaccuracies in the Kennicutt’s SFR 
determinations based on Her emission. 


6 THE ESCAPE ERACTION OF THE IONIZING 
PHOTONS 


In this Section, we investigate how many ionizing photons es¬ 
cape the galaxy volume. These photons co uld play a fundamen¬ 
tal role in the reion i zation of the Uni verse dVanzella et al. 1120121 : 
iKimm & Cen]l20l4IWise et al. 112014l) . 

Let the escape fraction of ionizing photons from sector i be 
defined as follows 




Qescii) 

Qstarsi^^ 


Fv{esc,i) dv 
^V\Ryd /'F 

Fv{stars,i) dv ’ 

Aisvrf 


( 2 ) 


where Qesed) denotes the rate at which ionizing photons leave sec¬ 
tor i through the outer edge (i.e. propagate above 2.3 kpc); Qsmrsii) 
is the rate of incident ionizing photons emitted by the stars con¬ 
tained in sector i; Fy(esc,i) and Fy(stars,i) are the correspond¬ 


Figure 13. Evolutionary dependences of the escape fraction of ionizing 
photons for sectors 1, 10 and 18 as well as for total escape fraction (see 
text). Also, the evolution of psc for the model M2 (Sector 18) is shown. 
The values of psc for the M2 model along Sectors 1 and 10 are » 0.01 at all 
ages, thus have not been shown. 


ing fluxes of ionizing radiation. It must be noted that Fy(stars) 
is not attenuated by the intervening gas; it is the overall flux pro¬ 
duced by all central stars. In practice, we perform an integration of 
both ionizing fluxes over the photon energy range (1 - 10)Ry for 
Fy(stars), and (1 - 180)Rv for Fy{esc), respectively. This is due 
to the fact that Fyistars, i) (calculated, as mentioned above, using 
Slarburst99) drops to zero beyond lORyd, but the escape radiation 
(Fy(esc, i)) contains also the photons at higher energies, emitted by 
the hot rarefied gas in the SWR. Notice also that high-energy Lyc 
photons are less absorbed, because the photoionization cross sec¬ 
tion decreases with photon energy increasing. 

The total escape fraction of the ionizing photons is defined as 

. _ QeAi) ... 

i.e. we simply sum up all the escape fractions, in all sectors. 

In Figure [T3] the evolution of for sectors 1, 10, and 18 as 
well as the total escape fraction is shown. As expected, lower val¬ 
ues of /ejc belong to low-number sectors, closer to the disk of the 
galaxy. The escape fraction of ionizing photons for the M2 model 
are in general very low (of the order of 0.01). Only the escape frac¬ 
tion for sector 18 (close to the symmetry axis) is significantly larger 
than 0, at least until 50 Myr. This is thus the only result we show 
for model M2 (see filled squares-fdot-dashed curve in Fig.ll3ll. 

We can notice from Figure[T3] that psc for Model Ml is on av¬ 
erage 0.4. This value is slightly hi gher than the value found in radia¬ 
tive h ydrodynamical simulations jKimm & Cen ll20l^ : IWise et al. I 
l2014h . On the other hand, model M2 produces much lower values 
for fgsc, lower than predicted by simulations and also lower than 
the values required to e nsure that DGs are significant for the reion¬ 
ization of the Universe. iFuiita et al, I ( 1200^ e.g. found that for DGs 
with SF efficiencies, i.e. the gas fraction converted into stars above 

fesc can be larger than 0.2 and this value is large enough for 
a significant contribution to the ionizing UV background. Notice 
however that the true SF efficiency in DGs could be smaller than 
6 %. In this case, DGs can not be the main responsible for the reion¬ 
ization of the universe. In the model L8F considered here, this SF 
efficiency is about 13%, thus our average of ~ 0.4 is consistent 
with Fujita’s calculations. 

As we have shown above, the small scale TDS, characteristic 
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of the M2 model, is necessary to reproduce the observed emission 
line spectra. Of course, in the presence of many TDSs in the DG, 
the total fesc would he further reduced. However, we shall not forget 
that our calculations are based on the outward-only approximation 
and it is possible that a fraction of the ionizing photons are emitted 
in non-radial directions. This would enhance the in model M2. 
It appears that some modifications in our modelling procedures (in 
particular, the relaxation of the outward-only approximation) are 
required in order to obtain values closer to realistic ones and 
less dependent on model assumptions. 

We also studied the changes of the energy distribution in the 
ionizing spectrum. For both Ml and M2 models, we show only re¬ 
sults relative to sector 1. In Figure [Tdb) we compare 1) the incident 
Lyc-spectrum from central stars in the energy range (13 - 140)eF 
for both models; 2) the spectrum of the radiation escaping the 
SWR; 3) the spectrum calculated at the outer edge of galaxy, pre¬ 
dicted by model Ml. The outer part of the galaxy absorbs very 
efficiently ionizing quanta in this photon energy range, because of 
the high density in Sector 1. Therefore, an outer ionization front is 
present in this sector. The Lyc spectra at the outer edge are very 
weak. Nevertheless, gaps in the Lyc spectra, mainly positioned at 
important ionization potentials of ions, are clearly distinguishable. 
Similar gaps were predicted by optimized photoionization mod- 
els of HII regions i n low-metallicity starburst galaxies (see e.g. 
iMelekh ||2006[|200% . As a result of these models, the optimal Lyc- 
spectra for HII regions in the blue compact galaxies SBS 0940-1-544 
and SBS 0335-052 was found. The gaps appear also in models of 
superwind bubbles jKoshmak & Meleldi1l2013l . r2014l) . 

As mentioned above, the incident Lyc spectrum from stars 
drops to zero at photon energies > l36eV (see Figures [T4h.b). but 
the emission of warm-hot gas contributes photons at higher ener¬ 
gies. The spectra in the range (0.01 - l)keV are shown in Figure 

It can be seen that the density distribution in the outer part 
of galaxy (from ’wall’ up to outer edge of the modelling volume) 
shapes the spectra up to photon energies of (0.4 - 0.5)keV. The 
spectrum at larger photon energies is practically undisturbed, be¬ 
cause of low the optical thickness for these photons. Therefore, the 
modelled spectrum in this energy range can be compared with X- 
ray observations. 

It must be also once again remarked that in our models we cal¬ 
culated the diffuse component of the the ionizing radiation in out¬ 
ward only approximation. Probably, the diffuse ionizing radiation 
does not play a crucial role in the formation of the global struc¬ 
ture of diffuse ionized gas. Nevertheless, as we have shown above, 
important emission lines can be partly formed in low-ionization 
zones of small spatia l scale. Existing 3D codes (such as Mocassin; 
lErcolano et al. II 2 OO 3 I 1 can calculate the diffuse ionizing field in re¬ 
gions shadowed by TDS, but they need a very large amount of com¬ 
puting time. Moreover, Mocassin is based on the statistical analysis 
of random propagation of photon packages. We think that for the 
problem at hand a deterministic approach, based on the detailed 
calculation of diffuse radiation propagation, is more useful. There¬ 
fore, we are developing an extension of the MPhM, able to follow 
in detail the diffuse ionizing radiation field (i.e. we will relax the 
outward-only approximation). Such method can be very useful for 
the correct investigation of the ionization structure in the shadow 
regions beyond the TDS. As mentioned above, this method can 
also help simulate the leakage of ionizing photons from a galaxy 
in a more realistic way. 
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a) Lyc spectra for photons energy (13-140)eV, age: 140 Myr. 
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b) Lyc spectra for photons energy (0.01 - 1)keV, age: 140 Myr. 
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Figure 14. Changes of the Lyc-spectrum from the stars during the transfer 
of ionizing radiation through the DG. The upper panel refers to the energy 
range between 13 and 140 eV. The lower panel extends up to 1 keV. All 
fluxes are given at a distance of 2.3 kpc from the center of galaxy (at the 
outer edge of the modelling volume). 


7 DISCUSSION 

In this paper we have described a way to combine detailed ChDS of 
DGs with a state-of-the-art MPhM, to study the emission properties 
of star-forming galaxies and compare them with observations. 

We have performed two types of analysis of the obtained 
MPhM results: (i) model predictions and (ii) consistency checks. 
Important emission line ratios have been calculated and compared 
with observations of a DG galaxy with structural parameters simi¬ 
lar to the simulated one (see Sect. I53t . The calculated escape frac¬ 
tion of ionizing photons fe,c is calculated, too, and compared with 
results of radiative transfer hydrodynamical simulations (Sect. nil. 
The results of these calculations are briefly summarized and dis¬ 
cussed in Sect. 17. II With consistency checks we mean instead tests 
in which some galactic properties are known because of the results 
of the ChDSs. The same galactic properties are re-derived as an 
observer would do, based on the MPhM results. The analysis of 
the oxygen abundances (Sect. lSAt and of the empirical determina¬ 
tions of the SFR in galaxies (Sect. [53] l belong to this category. The 
results of these two tests are briefly discussed below (Sect. IT^ . 
We shall not forget that our analysis is based on some simplifying 
hypotheses and there are potentially important physical processes 
(e.g. the dust reprocessing) that have not been considered. We will 
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briefly discuss these limitations in Sect. 17.31 Finally, some overall 
conclusions are drawn in Sect. [8] 


7.1 Model predictions 

We have calculated the intensities of observationally relevant emis¬ 
sion lines both integrating along angular sectors (Sect. and 
across apertures (Sect. Some line intensities and intensity ra¬ 
tios, such as [OIII]/15007/Flj8 and Htt/HyS, nicely fit into the range 
of observed intensity ratios (see Table 2) whereas some others do 
not. 

In particular, calculated [Oil] line intensities are very weak. It 
should be also remarked here that the average final oxygen abun¬ 
dance of the model L8F is 12-l-log(0/H)=7.84 (see RH13). We have 
also calculated model Ml at a higher metallicity, corresponding to 
the oxygen abundance 12-l-logO/H=8.2. For such model, the cal¬ 
culated [Oil] line intensities are still weak in comparison with ob¬ 
servations. This is due to the fact that the supershell created by the 
stellar feedback in ChDSs is not dense enough to block O-ionizing 
photons. O^ is thus very scarce (see Fig. m and this clearly implies 
very few O^ transitions. The inclusion of artificial density enhance¬ 
ments close to the supershell (model M2) helps creating a shell of 
singly ionized oxygen. For model M2, thus, also intensity ratios 
such as [OII]/13727///y3 are in agreement with observations. 

We will provide a comprehensive study on the models like Ml 
and M2 at different metallicities in a future paper. 

We have also calculated, for each angular sector, the escape 
fraction of ionizing photons This fraction is generally very 
high along directions close to the symmetry axis, i.e. towards 
the poles, because the superbubble propagates preferentially along 
these directions due to the very steep pressure gradients. The global 
fesc ranges from ~ 0.3 to ~ 0.6, with an average of about 0.4. 
Results of radiative hydrodynamical simulations of DGs generally 
show lower values of psc, although the results crucially depend on 
the details of the feedbac k scheme and on t he adopted star forma¬ 
tion efficiencies (see e.g. iFuiita et al. Il2003h . Escape fractions cal¬ 
culated for model M2 are instead generally low. They can be of the 
order of 0.3-0.5 only along directions close to the symmetry axis 
and only for a limited amount of time, otherwise their values are of 
the order of 0.01. With such low values of one can not expect 
that DGs signific antly contribute to the reionization of the universe 
as concluded by I Vanzella et al. I ( l2012h . In M2 models it was as¬ 
sumed that the TDS has the same density in all sectors. Of course 
this is a just a first-order approximation. TDS in different sectors 
would be characterized by different density values; in some sectors 
TDS might be even absent. Therefore, in reality, the Lyc photon 
escape should be higher in the vertical direction. Thus, evidently, 
more accurate calculations are necessary to address more in detail 
the issue of the escape fraction (see also Sect. [73] >. 

7.2 Consistency checks 

The employed ChDSs follow in detail the evolution in space and 
time of important chemical elements, in particular oxygen. Based 
on the chemical composition of different regions of the galaxy, the 
MPhM calculates intensity ratios, as extensively described in the 
previous sections. On the other hand, these intensity ratios can be 
measured in real galaxies and are the principal way observers can 
infer the chemical composition of galaxies. We study the results 
of our simulations as an observer would do and derive the oxy¬ 
gen abundances based on calculated intensity ratios. As explained 
above, we use the and R 23 diagnostic methods. 


In an ideal world, the abundances derived by means of inten¬ 
sity ratios should match the ’real’ ones, i.e. here the ones obtained 
from the ChDSs. This turns out to be the case only for some aper¬ 
tures and for some intervals of time (see Table 3). For other aper¬ 
tures and/or other intervals of time, the differences are significanf. 
These disagreements are due to two reasons: (f) the distribution of 
oxygen is very inhomogeneous and these inhomogeneities cannot 
be properly addressed by the applied diagnostic methods. These di¬ 
agnostic methods work well in relatively uniform regions, but fail 
in the presence of very inhomogeneously distributed chemical ele¬ 
ments. {ii) Diagnostic methods commonly assume that most of the 
oxygen is either singly or/and doubly ionized. This is not the case in 
model Ml, as we have shown in Fig. For this reason, the adop¬ 
tion of unresolved density enhancements in thin shells (model M2) 
helps in obtaining ’observed’ ionization stages of oxygen close to 
the ’real’ ones. 

In spite of that, our results indicate that the diagnostic methods 
should be used with a ’pinch of salt’ to say the less. For galaxies 
in which integral field specfroscopy is available, fhe difficulties of 
the diagnostic methods are reduced, as every chunk of the galaxy 
encloses a different abundance but the small-scale abundance inho- 
mogenei ties affect the integr ated spectral lines much less. Notice 
also that iLerov et al. I ( 1201 2h report of 1 kpc-scale variations and 
uncertainties of observed SFR indicators that are inherently e.g. 
caused by age effects and by intrinsic scatters of the indicators. 

The Htt emission from galaxies is extensively appli ed to de¬ 
rive their SFRs locally and globally (see lKennicutt |[T9^ and ref¬ 
erences therein), because it is in fact assumed to be the signpost 
of (short-living) massive stars. For an assumed IMF the number of 
massive stars is, therefore, correlated with the SFR. As we have 
shown (Fig. the total Ha emission is up to 60% larger than 
the one expect from the stellar ionizing flux of our model. If one 
takes into account that, in addition, ionizing photons are also lost 
by the leakage from the gala xy after a few tens of M yr (Fig. 1131 and 
that the Ha -SFR relation bv ICalzetti et al. I ( l2007h implies an even 
lower proportionality constant between SFR and Ha luminosity, 
one can conclude that the additional ionization by cooling photons 
from the warm-hot gas is almost equal to the stellar contribution. At 
later times, the Ha emission reduces, partly, due to the increasing 
amount of leaking ionizing photons. In this phase, an Ha -based de¬ 
termination of the SFR would underestimate the real galactic SFR 
(see again Fig. [Hi. However, given the uncertainties associated to 
this SFR indicator, a factor of two agreement can be considered as 
fair. 


Actually, it has been already reported in the literature that 
Ha could be a poor SFR indicator in DGs. In particular, it has 
been noticed that the Ha -based SFR estimates differ from fhe 
SFR estimates based on the UV emission at SFRs smaller that ~ 
0.01 Mq yr“’ and the dis crepancy can be larger than a factor of 
10 at SFR=10 Mq yr ' dLee et al. II2OO9I1 . This discrepancy has 
been associated to a top- light IMF in DGs with very mild SFRs 
dPflamm- Altenburg Il2009h . In such galaxies, there is some star for¬ 
mation but not enough generation of the very massive stars able to 
produce Ha emission. Notice, however, that our SFR (nominally 
2.67 • 10“^ Mq yr“') is not extremely low. In this SFR range the ef¬ 
fect of the top-light IMF should not be very significant. Notice also 
that the leakage of Ha photons can be significant for our model and 
it is surely even more significant for galaxies with lower masses, 
where galactic winds are more prominent (see RH13). It is thus 
reasonable to assume that t he discrepancy be tween Ha -based and 
UV-based SFRs noticed bv iLee et al. I d2009h also depends on psc- 
Notice however also that, as Lerov et al. “( boiTh point out. Star- 
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burst99 simulations of an evolving single stellar population imply 
an intrinsic scatter of almost 0.3 dex in Htr -based SFR estimates 
and about 0.5 dex in FUV-based determinations. We will add these 
issues in detail in a follow-up paper. 


7.3 Model limitations 

We have already outlined three limitations of the study presented 
in this paper, namely: 

• The MPhM is applied in a post-processing step, thus the en¬ 
ergy due to radiative processes does not feed back on the gas and 
does not affect the hydrodynamics. 

• The gas densities, obtained from the ChDSs, are always below 
a few particles per cm^, i.e. they never reach the peak densities of 
typical HII regions. 

• The MPhM calculations are based on the outward-only ap¬ 
proximation (the only one available with CLOUDY) 

We do not plan to relax the first simplification because this would 
require computationally expensive radiative hydrodynamics simu¬ 
lations. These calculations do not allow to explore a wide parameter 
space and, moreover, simplifying assumptions about the radiation 
field and the propagation of radiation need to be made. Notwith¬ 
standing the limitations, we think that our approach is reasonable 
because it is computationally inexpensive and the radiation trans¬ 
fer equations can be calculated with great accuracy. We do plan to 
relax the last two simplifications though. 

The second limitation can be relaxed by using a multi-phase 
ISM treatment, that takes into account (in the same computational 
cell) a warm-hot, diluted phase, and a colder, denser phase, that 
can reach the peak densities of typical HII regions and more. No- 
tice also that, for the galaxy NGC1569. IWestmoauette et al. I bOOTl 
I 2 OO 8 I) concluded from Ha line decompositions that the light profile 
is composed of a bright, narrow feature originating within a nar¬ 
row region of ~-700x500 pc in size, roughly centred on the bright 
super star cluster A, and of an underlying broad component emerg¬ 
ing from turbulent mixing layers on the surfaces the clouds which 
are embedded into the hot, fast-flowing winds from the young star 
clusters and experience evaporation and/or ablation of material. A 
new 3D chemo-dynamical code (cdFLASH) is already available 
(Mitchell et al., in preparation) and needs to be applied to the prob¬ 
lem at hand. 

The last limitation (the outward-only approximation) will be 
relaxed once a new CLOUDY wrapper (currently under develop¬ 
ment) will be available. 

It is, however, worth mentioning that the model presented 
here lacks another potentially relevant ingredient: dust. Dust can 
increase the opacities of many lines along all directions, thus af¬ 
fecting our results significantly. The reason why dust has not been 
included so far is not technical: the MPhM wrapper and the under¬ 
lying CLOUDY code can already treat dust. We have decided to 
neglect dust only for the purpose of clarity and in order to avoid 
cluttering of the (already long) text. We have however already per¬ 
formed a pilot study with dust. Assuming a gas-to-dust ratio of 
1400 (ten times larger than the Milky Way value, typical for DGs) 
we find very little differences between the results of a dusty model 
and the ones presented in this paper. However, for smaller gas-to- 
dust ratios the effect of dust can be very significant. We will explore 
the dust effect in a short follow-up paper. 


8 CONCLUSIONS 

The main conclusions of the present study can be summarized as 
follows: 

• Our combined ChDS-MPhM approach is able to produce in¬ 
tensity ratios of important emission lines which are in agreement 
with the observations. This is not the case for lines, which are 
severely underestimated. The reason is that the densities are on av¬ 
erage low and most of oxygen in the model galaxy is in high ion¬ 
ization stages. 

• An artificial density enhancement (simulating an unresolved 
thin shell) brings also [Oil] lines in agreement with observations. 

• We obtain typical escape fractions of ionizing photons (fesc) 
of the order of 40%. This escape fraction reduces to a few per cent 
if we consider the same density enhancements described above. 

• Diagnostic methods (T,, and R 23 ) are generally not able to re¬ 
produce the real abundance of heavy elements, obtained by means 
of the chemo-dynamical simulations. This is due in part to the in¬ 
herent limitations in the diagnostic methods, in part to the fact that 
in our simulation very little oxygen is neutral or singly ionized. 

• We have shown that the Kennicutt’s estimate of the SFR in 
galaxies based on the Ha emission is inaccurate. Part of the inac¬ 
curacy is due to the significant contribution of the warm-hot gas to 
the Ha emission. Part is instead due to the fact that the leakage of 
ionizing photons out of a typical starbursting galaxy can be signifi¬ 
cant and that naturally affects the global Ha emission. 
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